home *** CD-ROM | disk | FTP | other *** search
/ Nautilus 1992 July / Nautilus-3-8 / Nautilus-3-8.bin / Tools & Utilities / Techy Stuff / Source ƒ / sky source ƒ / MOON.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  10.7 KB  |  543 lines
  1. #include "sky.h"
  2.  
  3. /*
  4.  *    References:
  5.  *    Brown,
  6.  *    Improved Lunar Ephemeris
  7.  *    Supplement to A.E. 1968
  8.  *    Transformation corrections.
  9.  */
  10.  
  11. double k1, k2, k3, k4;
  12. double mnom, msun, noded, dmoon;
  13. double sinx();
  14. double cosx();
  15.  
  16. extern struct    moontab
  17. {
  18.     float    f;
  19.     char    c[4];
  20. } moontab[];
  21.  
  22. extern struct moont1
  23. {
  24.     float f1[2];
  25.     char c1[7];
  26. } moont1[];
  27.  
  28. moon()
  29. {
  30.     register struct moontab *mp;
  31.     register struct moont1 *mp1;
  32.     double dlong, lsun, psun;
  33.     double eccm, eccs, chp, cpe;
  34.     double q0, v0, t0, m0, j0, sn0, l0;
  35.     double arg1, arg2, arg3, arg4, arg5, arg6, arg7;
  36.     double arg8, arg9, arg10, arg11, arg12, arg13;
  37.     double arg14, arg15, arg16, arg17;
  38.     double arg18, arg19, arg20, arg21, arg22;
  39.     double dgamma, k5, k6;
  40.     double lterms, sterms, cterms, nterms, pterms, spterms;
  41.     double gamma1, gamma2, gamma3, arglat;
  42.     double xmp, ymp, zmp;
  43.     double temp1, temp2;
  44.     double shsd;
  45.     double obl2;
  46.     double planp[7];
  47.  
  48.     object = "moon        ";
  49.  
  50. /*
  51.  *    the fundamental elements - all referred to the epoch of
  52.  *    Jan 0.5, 1900 and to the mean equinox of date.
  53.  */
  54.  
  55.     dlong = 270.434164 + 13.1763965268*eday - .001133*capt2
  56.          + 1.9e-6*capt3;
  57.     argp = 334.329556 + .1114040803*eday - .010325*capt2
  58.          - 12.5e-6*capt3;
  59.     node = 259.183275 - .0529539222*eday + .002078*capt2
  60.          + 2.2e-6*capt3;
  61.     lsun = 279.696678 + .9856473354*eday + .000303*capt2;
  62.     psun = 281.220844 + .0000470684*eday + .000453*capt2
  63.          + 3.3e-6*capt3;
  64.  
  65.     dlong = fmod(dlong, 360.);
  66.     argp = fmod(argp, 360.);
  67.     node = fmod(node, 360.);
  68.     lsun = fmod(lsun, 360.);
  69.     psun = fmod(psun, 360.);
  70.  
  71.     eccm = 22639.550;
  72.     eccs = .01675104 - .00004180*capt;
  73.     incl = 18461.400;
  74.     cpe = 124.986;
  75.     chp = 3422.451;
  76.  
  77. /*
  78.  *    some subsidiary elements - they are all longitudes
  79.  *    and they are referred to the epoch 1/0.5 1900 and
  80.  *    to the fixed mean equinox of 1850.0.
  81.  */
  82.  
  83.     q0 = 177.481153 + 4.0923388020*eday;
  84.     v0 = 342.069128 + 1.6021304820*eday;
  85.     t0 =  98.998753 + 0.9856091138*eday;
  86.     m0 = 293.049675 + 0.5240329445*eday;
  87.     j0 = 237.352319 + 0.0830912295*eday;
  88.     sn0 = 265.869508 + 0.03345974279*eday;
  89.     l0 = 269.736239 +13.1763583100*eday;
  90.  
  91. /*
  92.  *    the following are periodic corrections to the
  93.  *    fundamental elements and constants.
  94.  *    arg4 is the "Great Venus Inequality".
  95.  */
  96.  
  97.     arg1 = 41.1 + 20.2*(capt+.5);
  98.     arg2 = dlong - argp + 33. + 3.*t0 - 10.*v0 - 2.6*(capt+.5);
  99.     arg3 = dlong + 3.*argp - 4.*lsun + 67. - 23.*t0 + 25.*m0;
  100.     arg4 = dlong - argp + 151.1 + 16.*t0 - 18.*v0 - (capt+.5);
  101.     arg5 = node;
  102.     arg6 = node + 276.2 - 2.3*(capt+.5);
  103.     arg7 = 152. + 119.*(capt+0.5);
  104.     arg8 = dlong + argp - 2.*lsun + 105. + t0 - 3.*q0;
  105.     arg9 = 313.9 + 13.*t0 - 8.*v0;
  106.     arg10 = dlong - argp + 112.0 + 29.*t0 - 26.*v0;
  107.     arg11 = dlong - argp + 21.*t0 - 21.*v0;
  108.     arg12 = dlong + argp - 2.*lsun + 273. + 21.*t0 - 20.*v0;
  109.     arg13 = dlong + argp - 2.*lsun + 303. + 8.*t0 - 12.*v0;
  110.     arg14 = 2.*lsun - 2.*node + 270. + 6.*t0 - 5.*v0;
  111.     arg15 = dlong + 2.*lsun - 3.*argp + 24.*t0 - 24.*v0;
  112.     arg16 = dlong - lsun + 262. + 12.*t0 - 15.*v0;
  113.     arg17 = dlong - lsun + 190. + 25.*t0 - 23.*v0;
  114.     arg18 = 43. - 8.*t0 + 15.*m0;
  115.     arg19 = node - lsun + 165. + 2.*m0;
  116.     arg20 = node + 290.1 - 0.9*(capt+.5);
  117.     arg21 = 115. + 38.5*(capt+.5);
  118.     arg22 = 216.0 - 8.*t0 + 15.*m0;
  119.     arg1 *= radian;
  120.     arg2 *= radian;
  121.     arg3 *= radian;
  122.     arg4 *= radian;
  123.     arg5 *= radian;
  124.     arg6 *= radian;
  125.     arg7 *= radian;
  126.     arg8 *= radian;
  127.     arg9 *= radian;
  128.     arg10 *= radian;
  129.     arg11 *= radian;
  130.     arg12 *= radian;
  131.     arg13 *= radian;
  132.     arg14 *= radian;
  133.     arg15 *= radian;
  134.     arg16 *= radian;
  135.     arg17 *= radian;
  136.     arg18 *= radian;
  137.     arg19 *= radian;
  138.     arg20 *= radian;
  139.     arg21 *= radian;
  140.     arg22 *= radian;
  141.  
  142.     dlong += (
  143.            0.84 *sin(arg1)
  144.          +  0.31 *sin(arg2)
  145.          +  0.04 *sin(arg3)
  146.          + 14.27 *sin(arg4)
  147.          +  7.261*sin(arg5)
  148.          +  0.282*sin(arg6)
  149.          +  0.04 *sin(arg7)
  150.          +  0.075*sin(arg8)
  151.          +  0.237*sin(arg9)
  152.          +  0.108*sin(arg10)
  153.          +  0.030*sin(arg11)
  154.          +  0.126*sin(arg12)
  155.          +  0.033*sin(arg13)
  156.          +  0.054*sin(arg14)
  157.          +  0.010*sin(arg15)
  158.          +  0.013*sin(arg16)
  159.          +  0.013*sin(arg17)
  160.          +  0.026*sin(arg18)
  161.          + 0.017*sin(arg19)
  162.          )/3600.;
  163.  
  164.     argp += (
  165.          - 2.10 *sin(arg1)
  166.          -  0.118*sin(arg4)
  167.          -  2.076*sin(arg5)
  168.          -  0.840*sin(arg6)
  169.           - 0.100*sin(arg7)
  170.          -  0.593*sin(arg9)
  171.          -  0.065*sin(arg18)
  172.         )/3600.;
  173.  
  174.     node += (
  175.            0.63*sin(arg1)
  176.          +  0.17*sin(arg4)
  177.          + 95.96*sin(arg5)
  178.          + 15.58*sin(arg6)
  179.          +  1.86*sin(arg20)
  180.          )/3600.;
  181.  
  182.     t0 += (
  183.         -6.40*sin(arg1)
  184.         -0.27*sin(arg7)
  185.         -1.89*sin(arg9)
  186.         +0.20*sin(arg22)
  187.          )/3600.;
  188.  
  189.     lsun += (
  190.         -6.40*sin(arg1)
  191.         -0.27*sin(arg7)
  192.         -1.89*sin(arg9)
  193.         +0.20*sin(arg22)
  194.         )/3600.;
  195.  
  196.     dgamma = -  4.318*cos(arg5)
  197.          -  0.698*cos(arg6)
  198.          -  0.083*cos(arg20);
  199.  
  200.     j0 +=
  201.            0.33*sin(arg21);
  202.  
  203.     sn0 +=
  204.          -  0.83*sin(arg21);
  205.  
  206. /*
  207.  *    the following factors account for the fact that the
  208.  *    eccentricity, solar eccentricity, inclination and
  209.  *    parallax used by Brown to make up his coefficients
  210.  *    are both wrong and out of date.  Brown did the same
  211.  *    thing in a different way.
  212.  */
  213.  
  214.     k1 = eccm/22639.500;
  215.     k2 = eccs/.01675104;
  216.     k3 = 1. + 2.708e-6 + .000108008*dgamma;
  217.     k4 = cpe/125.154;
  218.     k5 = chp/3422.700;
  219.  
  220. /*
  221.  *    the principal arguments that are used to compute
  222.  *    perturbations are the following differences of the
  223.  *    fundamental elements.
  224.  */
  225.  
  226.     mnom = dlong - argp;
  227.     msun = lsun - psun;
  228.     noded = dlong - node;
  229.     dmoon = dlong - lsun;
  230.  
  231. /*
  232.  *    solar terms in longitude
  233.  */
  234.  
  235.     lterms = 0.0;
  236.     mp = moontab;
  237.     for(;;) {
  238.         if(mp->f == 0.0){
  239.             mp++;
  240.             break;
  241.         }
  242.         lterms += sinx(mp->f,
  243.             mp->c[0], mp->c[1],
  244.             mp->c[2], mp->c[3], 0.0);
  245.         mp++;
  246.     }
  247.  
  248.     lterms +=
  249.         (294.e-9*eday - 9193./1296000.*dgamma + .0013)*sinx(1.0,0,0,0,2,0.)
  250.         +(-220.e-9*eday + 5282./1296000.*dgamma)*sinx(1.0,1,0,0,-2,0.);
  251.  
  252.     planp[1] = q0;
  253.     planp[2] = v0;
  254.     planp[3] = t0;
  255.     planp[4] = m0;
  256.     planp[5] = j0;
  257.     planp[6] = sn0;
  258.  
  259. /*
  260.  *    planetary terms in longitude
  261.  */
  262.  
  263.     mp1 = moont1;
  264.     for(;;){
  265.         if(mp1->f1[0] == 0.){
  266.             mp1++;
  267.             break;
  268.         }
  269.         lterms += sinx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
  270.             mp1->c1[2], mp1->c1[3],
  271.             mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
  272.         mp1++;
  273.     }
  274.  
  275.     lterms += sinx(-.189, 0,0,0,0, node) /*IAU*/
  276.         + sinx(-.013, -1,0,0,0, node) /*IAU*/
  277.         + sinx(-.013, 1,0,0,0, node); /*IAU*/
  278.     lterms += sinx(0.019, 0,0,0,0, 5.*t0-3.*v0+node+216.0);
  279.     lterms += sinx(0.016, 0,0,0,0, -5.*t0+3.*v0+node+075.0);
  280.     lterms += sinx(-.038, 0,0,0,0, 2.*node);
  281.  
  282. /*
  283.  *    solar terms in latitude
  284.  */
  285.  
  286.     sterms = 0.0;
  287.     for(;;) {
  288.         if(mp->f == 0.0){
  289.             mp++;
  290.             break;
  291.         }
  292.         sterms += sinx(mp->f,
  293.             mp->c[0], mp->c[1],
  294.             mp->c[2], mp->c[3], 0.0);
  295.         mp++;
  296.     }
  297.  
  298.     cterms = 0.0;
  299.     for(;;) {
  300.         if(mp->f == 0.0){
  301.             mp++;
  302.             break;
  303.         }
  304.         cterms += cosx(mp->f,
  305.             mp->c[0], mp->c[1],
  306.             mp->c[2], mp->c[3], 0.0);
  307.         mp++;
  308.     }
  309.  
  310.     nterms = 0.0;
  311.     for(;;) {
  312.         if(mp->f == 0.0){
  313.             mp++;
  314.             break;
  315.         }
  316.         nterms += sinx(mp->f,
  317.             mp->c[0], mp->c[1],
  318.             mp->c[2], mp->c[3], 0.0);
  319.         mp++;
  320.     }
  321.  
  322. /*
  323.  *    planetary terms in latitude
  324.  */
  325.  
  326.     pterms = 0.;
  327.     for(;;){
  328.         if(mp1->f1[0] == 0.){
  329.             mp1++;
  330.             break;
  331.         }
  332.         pterms += sinx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
  333.             mp1->c1[2], mp1->c1[3],
  334.             mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
  335.         mp1++;
  336.     }
  337.  
  338.     pterms +=
  339.           sinx(0.014, 0,0,0,0, -2.*t0+2.*v0+l0+285.0)
  340.          + sinx(0.027, 0,0,0,0, -1.*t0+1.*v0+l0+285.0)
  341.          + sinx(0.015, 0,0,0,0, t0-v0+l0+105.0)
  342.          + sinx(0.077, 0,0,0,0, 5.*t0-3.*v0+l0+215.6)
  343.          + sinx(0.025, 0,0,0,0, -6.*t0+4.*v0+l0+255.0)
  344.          + sinx(0.074, 0,0,0,0, -5.*t0+3.*v0+l0+051.6)
  345.          + sinx(0.018, 0,0,0,0, -4.*t0+2.*v0+l0+075.0)
  346.          + sinx(0.010, 0,0,0,0, -3.*t0+v0+l0+075.0)
  347.          + sinx(0.030, 0,0,0,0, 8.*t0-5.*v0+l0+125.0);
  348.     pterms +=
  349.          sinx(0.010, 0,0,0,0, -t0+2.*m0+l0+345.0)
  350.          + sinx(0.035, 0,0,0,0, 2.*j0+l0+168.0)
  351.          + sinx(0.018, 0,0,0,0, -2.*j0+l0+024.0)
  352.          + sinx(-.017, 0,0,0,0, l0)
  353.          + sinx(0.083, 0,0,1,0, 2.*node)
  354.          + sinx(0.215, 0,0,0,0, dlong) /*IAU*/
  355.          + sinx(-.012, -1,0,0,0, dlong) /*IAU*/
  356.          + sinx(0.011, 1,0,0,0, dlong); /*IAU*/
  357.  
  358. /*
  359.  *    solar terms in parallax
  360.  */
  361.  
  362.     spterms = 3422.700;
  363.     for(;;) {
  364.         if(mp->f == 0.0){
  365.             mp++;
  366.             break;
  367.         }
  368.         spterms += cosx(mp->f,
  369.             mp->c[0], mp->c[1],
  370.             mp->c[2], mp->c[3], 0.0);
  371.         mp++;
  372.     }
  373.  
  374. /*
  375.  *    planetary terms in parallax
  376.  */
  377.  
  378.     for(;;){
  379.         if(mp1->f1[0] == 0.){
  380.             mp1++;
  381.             break;
  382.         }
  383.         spterms += cosx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
  384.             mp1->c1[2], mp1->c1[3],
  385.             mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
  386.         mp1++;
  387.     }
  388.  
  389. /*
  390.  *    computation of longitude
  391.  */
  392.  
  393.     lambda = (dlong + lterms/3600.)*radian;
  394.  
  395. /*
  396.  *    computation of latitude
  397.  */
  398.  
  399.     arglat = (noded + sterms/3600.)*radian;
  400.     gamma1 = 18519.700 * k3;
  401.     gamma2 = -6.241 * k3*k3*k3;
  402.     gamma3 = 0.004 * k3*k3*k3*k3*k3;
  403.  
  404.     k6 = (gamma1 + cterms) / gamma1;
  405.  
  406.     beta = k6 * (gamma1*sin(arglat) + gamma2*sin(3.*arglat)
  407.          + gamma3*sin(5.*arglat) + nterms)
  408.          + pterms;
  409.     if(flags & OCCULT)
  410.         beta -= 0.6;
  411.     beta *= radsec;
  412.  
  413. /*
  414.  *    computation of parallax
  415.  */
  416.  
  417.     spterms = k5 * spterms *radsec;
  418.     hp = spterms + (spterms*spterms*spterms)/6.;
  419.  
  420.     rad = hp/radsec;
  421.     georad = 1.;
  422.     semi = .0799 + .272453*(hp/radsec);
  423.     if(dmoon < 0.)
  424.         dmoon += 360.;
  425.     mag = dmoon/360.;
  426.  
  427. /*
  428.  *    change to equatorial coordinates
  429.  */
  430.  
  431.     lambda += psi;
  432.     obl2 = obliq + eps;
  433.     xmp = georad*cos(lambda)*cos(beta);
  434.     ymp = georad*(sin(lambda)*cos(beta)*cos(obl2) - sin(obl2)*sin(beta));
  435.     zmp = georad*(sin(lambda)*cos(beta)*sin(obl2) + cos(obl2)*sin(beta));
  436.  
  437.     alpha = atan2(ymp, xmp);
  438.     delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
  439.  
  440. /*
  441.  *c    lunar eclipse computation
  442.  */
  443.  
  444.     shsd = 1.0183*hp/radsec - 969.85/rps;
  445.     temp1 = sin(shdecl)*sin(delta) + cos(shdecl)*cos(delta)
  446.      *cos(shra - alpha);
  447.     temp2 = atan2(sqrt(1.-temp1*temp1),temp1)/radsec;
  448.     temp2 = semi + shsd - temp2;
  449.     temp2 = temp2/(2.*semi);
  450.     if(temp2 >= 0.){
  451.         if(temp2 < 1.)
  452.             printf("Partial Lunar Eclipse: Mag. = %.4f\n", temp2);
  453.         else
  454.             printf("Total Lunar Eclipse: Mag. = %.4f\n", temp2);
  455.     }
  456.  
  457.  
  458.     geolam = lambda;
  459.     geobet = beta;
  460.  
  461.     geo();
  462.  
  463. /*
  464.  *    solar eclipse computation
  465. */
  466.  
  467.     if(!((flags&GEO)||(flags&HELIO))){
  468.         temp1 = sin(sundec)*sin(decl2) + cos(sundec)*cos(decl2)
  469.          *cos(sunra-ra);
  470.         temp1 = atan2(sqrt(1.-temp1*temp1),temp1)/radsec;
  471.         if(temp1 <= (semi2+sunsd)){
  472.             temp2 = (semi2+sunsd-temp1)/(2.*sunsd);
  473.             if(temp1 >= fabs(sunsd-semi2))
  474.                 printf("Partial Solar Eclipse: Mag. = %.4f\n", temp2);
  475.             else if(temp2 > 1.)
  476.                 printf("Total Solar Eclipse: Mag. = %.4f\n", temp2);
  477.             else
  478.                 printf("Annular Solar Eclipse: Mag. = %.4f\n", temp2);
  479.         }
  480.     }
  481.  
  482. /*
  483.  *    constants for occultation computations
  484.  */
  485.  
  486.     moonra = ra;
  487.     moonde = decl2;
  488.     moonsd = semi2;
  489.  
  490. }
  491.  
  492. double
  493. sinx(coef,i,j,k,m,angle)
  494. double coef, angle;
  495. {
  496.     double x;
  497.  
  498.     x = i*mnom + j*msun + k*noded + m*dmoon + angle;
  499.     x = coef*sin(x*radian);
  500.     if(i < 0)
  501.         i = -i;
  502.     for(; i>0; i--)
  503.         x *= k1;
  504.     if(j < 0)
  505.         j = -j;
  506.     for(; j>0; j--)
  507.         x *= k2;
  508.     if(k < 0)
  509.         k = -k;
  510.     for(; k>0; k--)
  511.         x *= k3;
  512.     if(m & 1)
  513.         x *= k4;
  514.  
  515.     return(x);
  516. }
  517.  
  518. double
  519. cosx(coef,i,j,k,m,angle)
  520. double coef, angle;
  521. {
  522.     double x;
  523.  
  524.     x = i*mnom + j*msun + k*noded + m*dmoon + angle;
  525.     x = coef*cos(x*radian);
  526.     if(i < 0)
  527.         i = -i;
  528.     for(; i>0; i--)
  529.         x *= k1;
  530.     if(j < 0)
  531.         j = -j;
  532.     for(; j>0; j--)
  533.         x *= k2;
  534.     if(k < 0)
  535.         k = -k;
  536.     for(; k>0; k--)
  537.         x *= k3;
  538.     if(m & 1)
  539.         x *= k4;
  540.  
  541.     return(x);
  542. }
  543.